File a12_regression diagnostics Topic Regression diagnostics See Class companion notes Data Abrasion Loss For a discussion of the data see the Visualizing Data by William Cleveland Due 3. One plot of your choice 4. The plot 5. The plot 7.2-7.5 One plot each 7.6 The last two plots 8.3 Two plots 9. Two plots Future: Additions: rotating 3d scatterplot Refinements: include units in labels 1. Get the abrasion loss data from the class web site. ##Run abrLossDat = read.table('abrasionLoss.txt') colnames(abrLossDat)=c('hardness','tenStrength','abrLoss') # reorder variables to put the dependent variable first abrLossDat = abrLossDat[,c(3,1,2)] abrLossDat # look at the values ##End The units of measure are Abrasion loss: g/hp-hour Tensile Strength: kg/cm**2 Hardness: degrees Shore 2. Define a function for looking at individual variables or regression residuals The smoothing parameter computed below has been multipled by 2 with the intent of oversmoothing ##Run eda1D = function(x,xlab=''){ oldpar = par(mfrow=c(2,2)) boxplot(x,ylab=xlab,main='Boxplot') qqnorm(x,ylab=paste(xlab,"Quantiles")) qqline(x) # stat = summary(x) # alternate smoothing parameter # smoothIQR = 1.5*(stat[5]-stat[2]) # could used and estimate smoothing parameter smoothPar = 2*bw.nrd(x) plot(density(x,width=smoothPar), xlab=xlab,ylab='Density', main='Density Plot', type='l') mtext(side=3,line=.3, paste('Smoothing parameter =',round(smoothPar,1))) hist(x,xlab=xlab) par(oldpar) return('done') } ##End 3. First 1-D distribution view of variables ##Run vnames = names(abrLossDat) for (i in 1:ncol(abrLossDat)){ windows() eda1D(abrLossDat[,i],xlab=vnames[i]) } ##End The windows, one per variable, are overplotted! Move the ones on top to see the ones below. In the QQplots note the suggestion of Abrasion loss: thick right tail Hardness: thin tails Tensile strength: thin tails 4. 2D Scatterplots ##Run library(lattice) windows() pairs(abrLossDat,panel=function(x,y){panel.smooth(x,y,pch=16,cex=1.2, col=3,lwd=2,)},gap=0) ##End In the top center panel y=Abrasion Loss and x= Hardness We see an almost linear relationship with low abrasion loss associaed with high hardness. In the top right panel y=Abrasion Loss and x= Tensile Strength We see an abrasion loss bump centered at about 160 for tensile strength. Mentally assess assessing what will happen with fitting about hardness and tensile strength using the two one variable plots is not so easy. The second row right panel y = hardness and x = tensile strength This panels shows the two variable domain for regression. Note the three points at the top left. Imagine thin rectangles to the right and below these points. There is an absence of data in such rectangles so there is little support for extending modeling results to these areas. Note there is little data in the lower left of the panel and that the two lowest hardness values are in the general location where the bump is located in the top right panel. This scatterplot matrix calls for brushing points with small hardness values to see where they are in all the panels. The Ggobi package could be installed to support brushing. The three points at the top left and the lowest hardness point with have a lot of leverage in the regression. That is, the regression may change substantially if all or some subset of these points were omitted. If all three of the values associated with each point are correct we will likely want to use them in the model. If there are some bad values, using such points can damage the model more than points in the center of the scatterplot. (the modeling domain). We can see high leverage point points pretty well in 2D domains but visualization become problematic in higher dimensional domains Scripts further below show calculation of point leverage. Classical regression assumes there is error in the dependent variable which is abrasion loss in this case. It assumes there is no error in the independent variables, which are hardness and tensile strength in this case. When there are errors in the independent variables statisticians called it the error in variables problem. Other disciplines encounter the same problem but use different names. There are no really good solutions to the problem, that don't involve new data. Some options are 1) Say the regression is conditional given the independent variable values 2) Use principle curve regression 3) Make some special assumptions. Apply singular value decomposition and use the left eigenvectors associated with largest eigenvalues as the independent variables. This hopes the contribution of bad values emerge primarily in the eigenvectors with low eigenvalues 5. A leverage example________________________________________ ##Run x = rnorm(50) y = 2*x+3 # Add high leverage points # and let y = -2x+3 for them x = c(x,-10,10) y = c(y,23,-17) #linear model # ~ read as: is modeled by # y is modeled by x model = lm(y~ x) summary(model) # model coefficients model$coef #(Intercept) x # 3.634926 -1.427085 plot(x,y,main="Linear regression with two bad high leverage points") abline(model$coef) ##End Two high leverage points have changed a perfect slope of 2 for the first 50 points into a slope of -1.43 6. An introduction to linear models and basic output In the linear modeling function lm() below data= indicates a data.frame with variables that can specified by name AbrLoss appear to the left of ~ which is read "is modeled by". AbrLoss is the dependent variable being modeled. The "." indicates that all the remaining variables in the the data.frame are to be used as predictor variables. By default the model also fits the grand mean. ##Run abrFit = lm(abrLoss ~ . ,data=abrLossDat) # summary of the fit abrSum = summary(abrFit) abrSum ##End Abrasion Loss Model Summary ___________________ Call: lm(formula = abrLoss ~ ., data = abrLossDat) Residuals: Min 1Q Median 3Q Max -79.385 -14.608 3.816 19.755 65.981 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 885.1611 61.7516 14.334 3.84e-14 *** hardness -6.5708 0.5832 -11.267 1.03e-11 *** tenStrength -1.3743 0.1943 -7.073 1.32e-07 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 36.49 on 27 degrees of freedom Multiple R-Squared: 0.8402, Adjusted R-squared: 0.8284 F-statistic: 71 on 2 and 27 DF, p-value: 1.767e-11 End Model Summary___________________________ Comment on the Model Summary Estimate column The fitted linear model is abrLoss = 855 -6.57*hardness-1.37*tenStrength Std. error column Estimate standard error t value Ratio of estimate to standard error Similar to a z-score with a theoretical mean of zero when the error degrees of freedom is 30 or large The Pr(>|t|) column indicate categories of statistical significance levels using a two sided t test. All three cooefficients have p < .001 Residual standard error: 36.49 on 27 degrees of freedom The residual standard error is 36.49. Perhaps other deterministic variables or stochastic components could drive this down. Multiple R-Squared: 0.8402 The Multiple R-Squared is ratio of the model predicted value sum of square about the grand mean to the data sum of square about the grand men. This can be computed from the residual sum of squares (RSS) from the two models. It is 1 - full_model_RSS/mean_only_RSS Modeling in the physical science often produces R-Squared values bigger than .5 (or 50% if converted to a percent) R-Squared values in the social sciences are often much less than .5. A high R-squared value simply indicates a reduction in sum of squares of residuals. It indicates neither causality, nor the applicability of the model to new data. One can often increase the R-squared value by generating a random variable and putting it in the model. For this reason some evaluation criteia include a penalty for the number linear predictors used in the model. F-statistic: 71 on 2 and 27 degrees of freedom, the p-value is 1.767e-011 The small p value for F statistics says the model is fitting the data a lot better than just the mean. This suggests that the model reduction in sum of squares is very unlikely to have happen at random assuming the predictor variables were unrelated to abrasion loss. Correlation of Coefficients: (Intercept) Hardness hardness -0.8335 tenStrength -0.7664 0.2992 The regression cooefficients are almost always correlated. Here Hardness and Ten.strength have a correlation of .2992 7. Common regression diagnostics See R. Dennis Cook and Sanford Weisberg "Regression Diagnostics With Dynamic Graphics", Technometrics 31(3) 1989 "Residuals and Influence in Regression", Chapman Hall 1982 and class notes. 7.1 Built in Diagostic Plots Today R and Splus has many built in diagnostic plots: ##Run !!! Click in the active plot to move to the next plot windows() plot(abrFit) ## Produces four plots (1) Residuals vs Fitted Values Candidate outliers are labeled (2) Standardized residual normal QQplot standardization divides residuals by their estimated variance asymmetry and thick tails are of concern (3) Spread location plot using y = sqrt(abs(standardized residuals)) x = fitted values (4) Standardized residuals versus leverage High leverage points strongly influence the fit If the point have bad y values the residuals still may be small. (Sometimes it is a bad x values the produce high leverage.) Sections below revisit some of these plots 7.2 Influential and high leverage cases________________________ ## Run abrInfluence = influence.measures(abrFit) infmat=abrInfluence$infmat influential=abrInfluence$is.inf nr = nrow(infmat) nc = ncol(infmat) rownam= rownames(infmat) colnam = colnames(infmat) windows(width=10,height=8) ngrps= (nr-1)%/% 4 + 1 pan = panelLayout(nrow=ngrps,ncol=nc+1, colSize=c(.7,4,4,4,4,4,4,4), rowSize=c(rep(4,ngrps-1),nr%%4), topMar=.5,bottomMar=.5) pan2 =panelLayout(nrow=1,ncol=nc+1, colSize=c(.7,4,4,4,4,4,4,4), topMar=.5,bottomMar=.5) for (i in 1:ngrps){ panelSelect(pan,i,1) subs = (4*i-3):min(4*i,nr) y =length(subs):1 panelScale(c(0,1),c(.3,length(y)+.7)) text(.5,y,rownam[subs],adj=.5,cex=.9) } for (j in 1:nc){ dat = infmat[,j] suspect = influential[,j] rx=range(dat) rx = mean(rx)+1.1*diff(rx)*c(-.5,.5) for (i in 1:ngrps){ panelSelect(pan,i,j+1) subs = (4*i-3):min(4*i,n) x = dat[subs] y = length(subs):1 col = ifelse(suspect[subs],"red","black") pch = ifelse(suspect[subs],16,1) cex = ifelse(suspect[subs],1.2,1) panelScale(rx,c(.3,max(y)+.7)) panelFill(ifelse(i%%2,'#D0D0D0','#C8C8FF')) points(x,y,col=col,pch=pch,cex=cex) panelOutline(col="white") if(i==ngrps)axis(side=1,tck=-.015,mgp=c(2,.3+.8*j%%2,0)) } } for (j in 1:nc){ panelSelect(pan2,1,j+1) panelScale() panelOutline() mtext(side=3,line=.3,colnam[j]) } pan3 = panelLayout(nrow=1,ncol=1, topMar=.5,bottomMar=.5) panelSelect(pan3,1,1) panelScale() title('Abrasion Loss Model: Influence Diagnostics') ## End The above layout could have vertical grid lines and vertical lines to indicated thresholds for taking a closer look at the data. Cases A and S have statistics flagged as filled red dots. 7.3 Large Studentized Residuals_____________________________ ##Run abrStRes = rstudent(abrFit) #studentized residuals x = seq(along=abrStRes) y = abrStRes fives = (x-1)%%5 + 1 plot(range(x),range(y),type='n', xlab='Case Sequence Number', ylab='Studentized Residuals', main="Unusual Residual Diagnostic Plot",las=1,) usr = par()$usr rect(usr[1],usr[3],usr[2],usr[4],col="#C0C0C0") abline(h=c(-2,2)) axis(side=1,col="white",tck=1,labels=F) box() col = c("red","orange","#00B000","#0080FF","purple") points(x,y,pch=16,col=col,cex=1.7) big = abs(y) > .9*2 text(x[big]+.02*diff(range(x)),y[big],as.character(x[big])) ##End 7.4 Spread-Location Plot_________________________________ Plot square root absolute residuals (y-axis) versus predicted values (x-axis) Look for horizontal smooth across the plot This may motivate transforming the data to stabilize the variance. In some data sets there in a nonlinear relationship between the the depandent variable and the predictor variables. A non-horizontal smooth can also be a symptom of nonlinearity. We expect some variation in the smooth and in takes some experience get comfortable about wiggles that are likley of little importance. There are more comments below. ##Run windows() scatter.smooth(abrFit$fitted,sqrt(abs(abrFit$res)), xlab='Fitted Abrasion Loss', ylab='Square Root Abs(Studentized Residuals)', main='Spread Location Plot',las=1,pch=16,cex=1.2,col="red") ##End A common regression concern is that variance of noise is not constant but rather increases (or decreases) with the magnitude of the dependent variable. The residual serve as a surrogate for noise and the predicted or fitted values serve as a surrogate for the underlying dependent variable values. The smooth in the plot does not look problematic. It is pretty flat. If there is a problem fairly often there is a data transformation to remove it. We don't have time pursue transforms but move on to other topics. 7.5 Normal QQ plot_______________________________________ Plot expected Normal Quantiles (x) versus sorted studentized residuals (y) Look for Tail thickness and extreme points *For a positive dependent variable examine plots for various values of Lamba in a Box-Cox transformation (not explained here) ##Run qqnorm(abrStRes,pch=16,col='red',cex=1.2) qqline(abrStRes) # "robust fit" ##End The left tail looks thick enough to cause some concern The right tail is a little thick. 7.6 Partial Residual or Added Variable Plots______________ Regular and Detrended Variables are usually included in models if their improvement to the fit is statistically significant. Some the improvement is to is due to one or a very few cases. This motivates examining such cases for the reliability of their values. Sometimes variables should be excluded from the models. Model deficiencies can also appear in the form of added variable outliers or high leverage points. The construction of an added variable plot is interesting. It regresses a selected predictor variable on all the other predictor variables to obtain a residual vector that is predictor variable's unique contribution to the full predictor variable vector space. The construction regresses dependent variable on all the other variables to obtain a residual vector that represents what has not been explained by the other predictor variables. The last construction step regresses the dependent variable residual vector on the predictor variable residual vector. The slope cofficient for this model is the same as the coefficient the predictor variable when fitting all the predictor variables. However the points in a scatterplot of y = dependent variable residual vector versus x = predictor variable residual vector reveals which cases are important to fitting what has not yet been modeled. The procedure is typically applied all individual predictor variables when there are not too many. Some of the linear regression can be avoided as the example script illustrates. Since it is easier to visaully assess departures from a horizontal line than from a regression line line with non-zero slope, the script also also produces detrended added variable plots that remove the regression fit. Adding smoothes to this plot helps us seen a possible function relationships in the residuals. ##Run # Create a matrix that omits the dependent variable depname=colnames(abrLossDat)[1] mat = as.matrix(abrLossDat[,-1]) vnames = colnames(mat) for (i in 1:ncol(mat)){ windows() xres =lm(mat[,i]~mat[,-i])$res scatter.smooth(xres,abrFit$res, xlab=paste('Adjusted',vnames[i]), ylab='Abrasion Loss Residuals', main='Detrended Partial Residual Plot', col="red",pch=16,cex=1.2) } for(i in 1:ncol(mat)){ windows() plot(xres,abrFit$res+abrFit$coef[i+1]*xres, xlab=paste('Adjusted',vnames[i]), ylab=depname, main='Partial Residual (Added Variable) Plot', col="red",pch=16,cex=1.2) abline(0,abrFit$coef[i+1]) } ##End There are four overplotted windows. Move the top ones to see the ones below Both added variable plots show many points supporting the regression lines. This is good. The quadratic appearance in the partial residual variable for tensile strength plots is suggestive but what action to take is not clear. 8. Surface and Contour Views of a Loess Smooth ## Run abrSmooth = loess(abrLoss~hardness*tenStrength, degree=2,span=.5,data=abrLossDat) # generate 2D grid hardness=abrLossDat$hardness tenStrength=abrLossDat$tenStrength gridMargins = list( hardness=seq(min(hardness),max(hardness),length=20), tenStrength=seq(min(tenStrength),max(tenStrength),length=20)) gridFull = expand.grid(gridMargins) # domain for prediction # Generate predict values for grid values # Note that extrapolation may be involved in some corners. # Ideally the plotting is restricted to the convex hull # based on the available data. abrGridPredict = predict(abrSmooth,gridFull) ## End 8.1 Contour Plot ##Run contour(gridMargins$hardness,gridMargins$tenStrength,abrGridPredict, nlevels=15, xlab='Hardness',ylab='Tensile Strength',main='Predicted Abrasion Loss') ##End You many need to close or move windows to bring the active window into view. 8.2 Perspective View ##Run persp(gridMargins$hardness,gridMargins$tenStrength,abrGridPredict, xlab='Hardness',ylab='Tensile Strength',zlab='Abrasion Loss') title(main='Predicted Surface - Try 1') ##End 8.3 Lattice (Trellis) Graphics for surfaces ##Run gridDat = data.frame(Abr = as.vector(abrGridPredict), Hardness=gridFull$hardness, TenStr=gridFull$tenStrength) library(lattice) # The followingsequence of 5 windows rotates the 3D plot # Move the windows to see the progression. # # One student animated the 3D plots in R but perhaps used # a different 3d plotting functions that wireframe. # I have not gotten wireframe to run in a for() loop. windows() wireframe(Abr~Hardness*TenStr,data=gridDat,drape=T) windows() rotZX = list(z=55,x=-60) wireframe(Abr~Hardness*TenStr,screen=rotZX,data=gridDat,drape=T) windows() rotZX = list(z=70,x=-60) wireframe(Abr~Hardness*TenStr,screen=rotZX,data=gridDat,drape=T) windows() rotZX = list(z=105,x=-60) wireframe(Abr~Hardness*TenStr,screen=rotZX,data=gridDat,drape=T) ## End Now a level and a contour plot windows() levelplot(Abr~Hardness*TenStr,data=gridDat,main="Lattice Level Plot") # Nice colors could be added as in previous a previous assignment windows() contourplot(Abr~Hardness*TenStr,data=gridDat,cuts=15,main="Lattice Contour Plot") ##End The juxtaposition of a contour plot belows a surface plot can be helpful. 9. Conditioned Plots for Surface Examination Conditioning partitioning data into subsets based on value of categorical variables and/or intervals of continuous variables. The enables two tasks 1) Focused assessment of the individual subsets. The assessment can include models or smooths and residuals. It often substantially the restricts the variation due to the conditioning variables It lowers the plotting dimensionality 2) Provides a basis for organized comparison across subsets There are many 3 variable applications for which contour and surface plots are sufficient for the task at hand. Conditioning is occational helpful. Conditioning becomes more important with there are four or more variables. The conditioned maps had five variables: two geospatial coordinates, and three attributes with one show in color two addressed by conditioning. The example coplot() example below only illustrates conditioning on one variable but coplot() will handle two conditioning variables. to The co.intervals() functions below define intervals for continous variables. The interval are allowed to overlap. If so they do not provide a strick partitioning. In the conditioned choropleth maps I choice to use strict partitioning to simplify explanation and to enable analysis of variance like comparisons across subsets. Cleveland using overlapping intervals, sometime refered to as shingles, in a smoothing context. Loess smooths are less reliable at the edges so having a few more points at the edges is helpful. ## Run hardnessIntervals = co.intervals(hardness,number=4,overlap=1/4) tenStrIntervals = co.intervals(tenStrength,number=4,overlap=1/4) coplot(abrLoss~hardness | tenStrength, data=abrLossDat, given.values=tenStrIntervals, panel=function(x,y,...) panel.smooth(x,y,span=.7)) ## End The top panel shows the conditioning intervals for tensile strength. The two lowest bars (position and value wise) give the conditioning intervals for the lowest panels in the matrix of panels below. The increasing of intervals corresponds to reading the matrix of scatterplot left to right, starting with the bottom row and working up. At a first glance the conditioned smooth look linear and differs mainly in the intercept. To appreciate the changes in intercept compare the vertical distances to the red lines from a common grid location in each panel such as x=50, y=300. In the bottom left panel the red line is a little above the point. In the top right panel it is far below the point. The last three points in the lower left plot suggest a change in slope. ## Run windows() coplot(abrLoss~tenStrength | hardness, data=abrLossDat, given.values=hardnessIntervals, panel=function(x,y,...) panel.smooth(x,y,span=.8,col="black")) ##End The smoothes against tensile strength conditioned on hardness show more variety. Note that the conditioned plot shows data points and smooth so suggest where residual are large, at least conditionally. Look at a smoothed surface plot again. It does not show the data. A surface plots plus points with line from points to the surface can show large residuals. A translucent surface can help reveal points and lines above and below the surface.